Biasing the center of charge in molecular dynamics simulations with empirical valence 
bond models: free energetics of an excess proton in a water droplet 
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Multistate empirical valence bond (EVB) models provide an accurate description of the energetics 
of proton transfer and solvation in complex molecular systems and can be efficiently used in molecular 
dynamics computer simulations. Within such models, the location of the moving protonic charge can 
be specified by the so called center of charge, defined as a weighted average over the diabatic states 
of the EVB model. In this paper, we use first order perturbation theory to calculate the molecular 
forces that arise if a bias potential is applied to the center of charge. Such bias potentials are 
often necessary when molecular dynamics simulations are used to determine free energies related to 
proton transfer and not all relevant proton positions are sampled with sufficient frequency during the 
available computing time. The force expressions we derive are easy to evaluate and do not create any 
significant computational cost compared with unbiased EVB-simulations. As an illustration of the 
method, we study proton transfer in a small liquid water droplet consisting of 128 water molecules 
plus an excess proton. Contrary to predictions of continuum electrostatics but in agreement with 
previous computer simulations of similar systems, we observe that the excess proton is predominantly 
located at the surface of the droplet. Using the formalism developed in this paper, we calculate the 
reversible work required to carry the protonic charge from the droplet surface to its core finding a 
value of roughly 4 ksT. 



I. INTRODUCTION 



Proton transfer is of crucial importance for a variety 
of processes in nature and technology ranging from ATP 
synthesis in living cells [l[ and enzymatic catalysis [^ 
to chlorine chemistry on stratospheric ice particles in- 
volved in polar ozone depletion [1, 0, Q and electrical 
power generation in hydrogen fuel cells [g, 0]. As pro- 
ton transfer involves cleavage and formation of covalent 
bonds its computer simulation is challenging. Ab initio 
methods such as density functional theory or wavefunc- 
tion based methods can model chemical reactivity but 
are computationally extremely expensive [^, Q . Molec- 
ular dynamics trajectories obtained with such methods 
are short and the number of collected transfer events is 
not sufficient to carry out a thorough statistical analysis 
capable of revealing the details of the mechanism. Re- 
cently, however, computationally efficient empirical va- 
lence bond (EVB) models for proton transport have been 



proton transfer events occur. 

Empirical valence bond models have been used to 
study proton solvation and transport in many systems 

13, water clusters [21,22,2^, 

[2y, [231, acidic aqueous solu- 

the water- vapor interface [28|, and biolog- 



including bulk water 
water filled pores 
tions 113, Hr 






developed by Voth and collaborators 




based on pio- 
These models. 



and Vuilleumier and Borgis 
neering work of Warshel [1^ |19, 
in which the Born-Oppenheimer surface is obtained as 
the lowest root of a secular equation involving empiri- 
cally modeled diabatic states and coupling elements, ac- 
curately describe the energetics of bond cleavage and for- 
mation but are computationally far less expensive than 
ab initio methodologies. With forces computed with the 
Hellmann-Feynman theorem one can perform nanosec- 
ond molecular dynamics simulations during which many 



ical systems [18, 19, 20, 29, 30i]- Many of these appli- 
cations have been surve yed in a recent review article by 
Voth and collaborators [ij|. A quantitative understand- 
ing of proton transport and solvation often requires the 
calculation of free energies related to different positions 
of the excess proton. For instance, the rate for proton 
translocation across membrane nanopores is mainly de- 
termined by the free energetic cost required to remove 
the proton from the bulk and bring it into the interior 
of the pore [2^. The definition of the position of the 
excess proton, however, is not straightforward in empir- 
ical valence bond models. In bulk water, for instance, 
the excess proton occurs in a variety of different con- 
figurations with the so called Eigen and Zundel cations 
as the limiting cases [3, [ill, [31| ■ Whereas in the Eigen 
cation, (11904)+, the proton is strongly associated with 
one particular water molecule that donates three con- 
tracted hydrogen bonds to adjacent water molecules, the 
excess charge is equally shared by two water molecules 
in the Zundel cation, (H502)+. Under the influence of 
the fluctuating hydrogen bond network, these structures 
freely convert into each other by proton hops along hy- 
drogen bonds on a picosecond time scale. Due to this 
fluxional character of the hydrated proton, identification 
of a particular proton as the excess charge is ambiguous. 



In the empirical valence bond model this ambiguity 
can be resolved by defining the center of charge, the 
weighted average of the positions of the excess charge 
in the diabatic states [32] . The center of charge provides 
a meaningful definition of the proton position in a way 
that changes continuously in time. The free energetics of 
proton translocation can then be determined in molec- 
ular dynamics simulations by accumulating statistics on 
the position of the center of charge. This straightforward 
approach works well if the free energy differences do not 
exceed the thermal energy k^T appreciably. If they do, 
the molecular dynamics simulation most likely fails to ex- 
plore all important regions of configuration space within 
the available computing time. In such cases, it is pos- 
sible to enhance the sampling of configuration space by 
introducing appropriate bias potentials that guide the 
system towards configurations that would not be sam- 
pled otherwise |[33|. For instance, in the case of proton 
transfer through a membrane, penetration of the proton 
into the pore interior is observed in a molecular dynamics 
simulation only if the related desolvation penalty is com- 
pensated by a suitable bias. Correction for the bias then 
permits to deduce ensemble averages and free energies 
for the bias-free system. 



When bias potentials are used in molecular dynamics 
simulations it is necessary to calculate the atomic forces 
resulting from the bias. In the case of a bias potential 
acting on the center of charge, such a calculation requires 
to determine the derivatives of the center of charge with 
respect to the particle coordinates. How to do that using 
first order perturbation theory is the central subject of 
this paper. The expressions obtained in this way permit a 
computationally inexpensive evaluation of the bias forces 
because they depend only on quantities already required 
in the calculation of regular EVB-forces. As an example, 
we apply the formalism developed in this way to study 
the solvation of an excess proton in a small water droplet. 
We note that free energies as a function of the position 
of the center of charge have been calculated before in bi- 
ased EVB-simulations [27l. [sil. [35| . In these publications, 
however, the calculation of the forces resulting from the 
bias is not explained and it is unclear how it is done. 



The remainder of this article is organized as follows. 
In the next section we will briefly review the multi-state 
empirical valence bond model and define the center of 
charge. First order perturbation theory is then used in 
Sec. mil to derive the atomic forces resulting from a bias 
potential acting on the center of charge. To illustrate the 
formalism, we calculate the reversible work required to 
transfer an excess proton from the surface to the interior 
of a small water droplet defined in Sec. IIVI and present 
the results in|Vl Some conclusions are given in Sec. IVII 



II. MULTI-STATE EMPIRICAL VALENCE 
BOND MODEL 

In this section we briefly review the multi-state em- 
pirical valence bond model in order to introduce the ter- 
minology and set the notation. For a more detailed de- 
scription we refer the reader to the original publications 
[lOl [l5| . The system, which includes of a certain num- 
ber of water molecules plus an excess proton and pos- 
sibly some other species, is described by 3A^ atomic co- 
ordinates X = {x^ ,x^, ■ ■ ■ ,x^^}. The position of par- 
ticle i is also denoted by r.^ such that we can write 
X = {ri,r2, • ■ • jTtv} when it is convenient. In the em- 
pirical valence bond model one imagines that the elec- 
tronic state \ip) for a particular atomic configuration is a 
superposition of L diabatic valence bond states \ipi), 



W ^^C^\ipi), 



(1) 



where the Ci are real expansion coefficients. In each state 
chemical bonds are assigned such that a different oxygen 
atom holds the hydronium ion HgO"*". In some states 
this assignment may correspond to a rather distorted ge- 
ometry of the hydronium ion and the neighboring water 
molecules. Such states will have a high energy and there- 
fore contribute only little to the ground state. Only a 
small number of states are considered which are chosen 
to account for different plausible routes for proton trans- 
fer. To construct these states, one starts from a pivot 
state and then includes other states that are accessible 
[l(j |. Typically, about 10 states are required to accurately 
model proton transfer in bulk water. 

To find the Born-Oppenheimer energy surface as a 
function of the atomic coordinates one does not de- 
termine the electronic state \ip) from first principles. 
Rather, all L^ matrix elements Hij{x) = {(pi\H\ifi) of 
the electronic Hamiltonian H are modeled empirically as 
functions of the atomic coordinates x. While the L di- 
agonal elements are the potential energies corresponding 
to particular hydronium ions with fixed chemical bond- 
ing topology, the L{L — l)/2 non-diagonal elements pro- 
vide the coupling between the diabatic states that enables 
chemical reactivity. Postulating that the diabatic states 
are orthogonal to each other, solution of the eigenvalue 
problem 



Hci = Eid 



(2) 



yields the adiabatic ground state. Here, c^ = 
{cjo, Cii, • • • , Ci(2,_i)} is the eigenvector of dimension L 
belonging to the eigenvalue Ei and its elements are the 
corresponding coefficients in the superposition of Equ. 
^ . (Here and in the following we number states starting 
with 0.) The Born-Oppenheimer potential energy surface 
is then given as the expectation value of the energy in the 
ground state. 



Eo{x) 






t,j 



coiCQjH,j{x), 



(3) 



where coi is the coefficient of state i in the ground state. 
The corresponding forces can be calculated using the 
Hellmann-Feynman theorem [3a |: 



coefficients Coi depend on the particle coordinates we ob- 
tain 






dx°' 



(4) 



dqP 
dx' 



2 2^ ^0^9^'?^ +2^ Co, 



7J37 = 2> co,Tj— '7r + > Co,77--. (8) 






./3 



Here, F" is the force acting on degree of freedom a. This 
procedure, which retains quantum mechanics on a rudi- 
mentary level, effectively interpolates in a smart way be- 
tween the diabatic states for which empirical potential 
energy surfaces for fixed chemical bonding topology are 
known. 

We can now define the center of charge Qc as the 
weighted average of the hydronium oxygen position over 
all diabatic states. 



yi cg,q,; 



(5) 



Here, q^ is coordinate (3 of the hydronium oxygen of di- 
abatic state i. The partial derivative dq[ /dx'^ is unity 
only if x" is the coordinate [3 of the oxygen carrying the 
hydronium in state i and it vanishes otherwise. 

According to Equ. ([8]), we need to evaluate deriva- 
tives of the elements of the ground state eigenvector cq 
with respect to the coordinates of all particles. This can 
be accomplished with ordinary first order perturbation 
theory (see, for instance, Ref. [33|). To do so, we first 
consider how the matrix element Hij of the Hamiltonian 
changes if one particular coordinate x" is changed by an 
infinitesimal amount A, i.e., 



Hi 



Hij + A— 



Here, q^ is the position of the hydronium oxygen in dia- 
batic state i and €% is the statistical weight of the dia- , , ■ ^ j j 
'^\. , . , „,.,„.. where we have introduced 



= Hi 



AA", 



batic state i in the adiabatic ground state. This definition 
of the center of charge takes into account the delocalized 
nature of the protonic charge and can be used to follow its 
migration along the fluctuating hydrogen bond network. 



III. BIAS FORCES 

We now imagine that a bias potential VeCqc), which 
depends on the center of charge qc, is added to the po- 
tential energy of the system. (If the bias potential is a 
function of other degrees of freedom in addition to qc, the 
formalism below needs to be adapted appropriately.) De- 
pending on the particular situation one wants to study, 
the bias potential VB(qc) is designed to control the po- 
sition of the protonic charge and enhance the sampling 
of the configuration space regions of interest. If one uses 
such a function in a molecular dynamics simulation it is 
necessary to determine the corresponding forces via 



D" ^ TT 

'=- dx^ ''' 



(9) 



(10) 



Note that the derivatives Df- are already available from 
the regular EVB force calculation. 

To calculate the derivative of coi with respect to x°' we 
need to determine how coi changes if we add the pertur- 
bation AD" to the Hamiltonian H . (Here, D" denotes 
the matrix consisting of the elements Df,.) Then, we can 
calculate the derivative from: 



dcpi coi(A)-coi 

- — = hm ; 

dx'^ A^o A 



(11) 



F^{^) = 



dVBinc 
dx°' 



(6) 



The argument in coi(A) indicates the dependence of coi on 
the strength of the perturbation, i.e., on the displacement 
of coordinate a;". If no argument is written for coi the 
value for a vanishing perturbation strength is implied. 

The ground state coefficients Coi(A) as a function of 
the displacement A can be calculated by expanding the 
coefficients and the energies in powers of A and truncating 
the expansion after first order. 



In the following we will derive an explicit expression that 
permits to calculate exactly these forces. 

Using the chain rule and the definition of the center 
of charge from Equ. ([5] )we can write the derivative of 
Vb{cIc) with respect to the coordinate a;" as 



Ek{X) 



CkiW 




0{X% 



(12) 



3^k 



dVBiqc) 
dx" 



3 

E 

/3=1 



dVBicic) 



7^ dx'- ' 



(7) 



where the superscript denotes the order of the term and 
A{X) is a normalization constant. Standard time inde- 
pendent perturbation theory then yields [38| 



where gf is coordinate f3 of the center of charge qc. To 
calculate the derivatives dq^/dx°' of the center or charge 
with respect to the particle coordinate we differentiate 
Equ. (O. Taking into account that also the ground state 



XE^. - y CklXD",^Ckm, 



A«(i) - 



E 



(0) 



E 



(0) 



(14) 
(15) 



Since the normalization factor A{X) is unity to first order, 
we obtain 



XEkiX) - Ei°^+\Y,Dr,nCkiCkm, 

Lrrt 
J^k ^k ~^j 

Inserting Equ. p7|l into Equ. (Ilip we obtain 



dco 






dx' 



j^O l,m 



En — E-i 



i-ji, 



(16) 
CH- (17) 



(18) 



where, for simphcity, we have omitted the superscripts 
for the unperturbed energies. Thus, we now have all 
elements in place necessary to calculate the atomic forces 
resulting from the bias potential using Equs. ([7]) and ([5|). 
When evaluating the bias forces according to Equ. ([7]) 
in a simulation it is important to note that derivatives 
dcoi/dx" of the ground state coefficients appear only in 
contraction with the vector coiqf (see Equ. ([8])). Insert- 
ing Equ. (fT8|) into Equ. ([8]) we obtain 



9x" 



2r"^ + V c^ ^ 



(19) 



where for convenience T"^ is defined as the multiple sum 



na/3 






En — Ei 



(20) 



The quadruple summation in the above equation is the 
compuationally most expensive part of the calculation of 
the bias forces (compared to the overall cost of an EVB 
simulation, however, the cost of the bias force calculation 
is negligible). 

The number of operations required for the calculation 
of r"^ now depends on the order in which the summa- 
tions are carried out. This situation is similar to that of- 
ten encountered in electronic structure calculations (see, 
for instance, Ref. [33|). One may, for instance, decide to 
carry out the summation over the indices m and I first 
for all j , 



B" =2j'^'™^0mCj7, 



and do summation over j for all i after that. 



(21) 



dx 



j¥o 



Eq — Ej 



(22) 



Finally, V^^ is calculated by summing over i. 



dcn 






dx" 



(23) 



The total number of operations of this procedure is of 
order L^ (recall that L is the number of diabatic states). 
While most summation orders lead to the same L^- 
scaling, there is a particular way to carry out the summa- 
tions for which the number of required operations scales 
as L^. This is achieved by first summing over i for all j. 



Qf = E 



CjiCoiQ'^ 



P 



Then, one carries out a summation over j for all /, 






CjiQj 
Eq - Ej 



The last step consists of a summation over m and 



l.tn 



(24) 



(25) 



(26) 



In this case, each of the three steps requires of the order of 
L^ operations so that also the total number of operations 
is of that order. 

Using these expressions one can calculate the forces 
caused by the bias potential at little extra cost and use 
them in a molecular dynamics simulation. Due to the 
denominator containing the energy difference between 
the ground state energy and higher eigen-energies of the 
EVB-Hamiltonian the above expressions are valid only in 
the non-degenerate case. Although the degenerate case 
can be treated with slightly modified expressions, this 
has never been necessary in our simulations. 

We close this section on the formalism by mentioning 
an alternative viewpoint that, however, yields exactly the 
same expressions for the bias force. This perspective is 
based on the observation that the center of charge can be 
expressed as a derivative of the (appropriately modified) 
total energy with respect to the components of constant 
electric field £ = {£^ , S^ , E^} which couples to the hydro- 
nium oxygens of the diabatic states. More specifically, 
the elements of the Hamiltonian matrix are modified in 
the following way: 



Hij {x, £) = Hij (x) + Sij {£ ■ qi 



(27) 



Then, the Hellmann-Feynman theorem can be used to 
show that the components of the center of charge are 
given by derivatives with respect to the field. 



dEo{x,£) 
d£P ■ 



(28) 



Using this relation, the derivatives of the center of charge 
with respect to the particle coordinates can be expressed 
as 



9gf _ d^Eo{x,£) 
dx<^ dx<^d£P 



d^Eo{x,£) 
d£Pdx<^ 



"d£P' 



(29) 



where we have exploited that the differentiation order 
can be exchanged and the derivative with respect to £^ 



is evaluated at f = 0. The derivatives of the forces F" 
with respect to the particle coordinates required in the 
above equation can be computed from Equ. ^ using 
perturbation theory. In contrast to the development de- 
scribed above, the perturbation is done with respect to 
the 3 electric field components rather than to the 3A^ 
particle coordinates. Nevertheless, this approach results 
in exactly the expressions of Equs. ^, P^ . and (PO)) . 
which can be evaluated efficiently as explained above. 



TABLE I: 


Parameters for the 


runs 


with biasing 


potential. 


fcbias 




Ry 


)ias 


length 



10.0 kcal/mol A" 
6.0 kcal/mol A" 
6.0 kcal/mol A' 
3.0 kcal/mol A" 



0.0 A 

1.5 A 

3.0 A 

4.0 A 



5 ns 
5 ns 
5 ns 
5 ns 



IV. SYSTEM AND SIMULATIONS 

As an illustration of the method described in the pre- 
vious section, we apply it to a system consisting of a 
cluster of Nw = 128 water molecules plus one excess 
proton. Thus, the total number of atoms is A^ = 385. 
The interaction energies and forces were calculated us- 
ing the multi-state empirical valence bond model of Voth 
and collaborators [lOl, llli, which was demonstrated to 
accurately describe proton transfer in bulk water. The 
simulations were carried at a temperature of T = 280 
K controlled using an Andersen thermostat [39| in which 
particle momenta are randomly assigned from a Maxwell- 
Boltzmann distribution with a rate of i^ = 0.1 ps~^. The 
equations of motion were integrated using the velocity 
Verlet algorithm with a time step of At — 0.5 fs for a 
hydrogen mass of 2 amu and an oxygen mass of 16 amu 
(a hydrogen mass of 2 rather than 1 amu leaves the struc- 
tural properties of the system unchanged, but allows for 
a larger time step). At each time step, the diabatic states 
are determined using the algorithm developed by Schmitt 
and Voth [I^, [ill such that the set of states is contin- 
uously adapted according to the location of the excess 
charge. The EVB-Hamiltonian is diagonalized using the 
Jacobi method for symmetric matrices. Typically, the 
number of diabatic states of the empirical valence bond 
model fluctuates between 8 and 10. 

To prevent single water molecules to evaporate from 
the cluster we have enclosed the whole system in a con- 
fining potential of the form 



VconfCa;) 



fc, 



N 
conf V^/^cm 



J2{rr - RcontfOirr - Rcont), (30) 



where 9{x) is the Heaviside step function and 



I "^i -tvcni I 



(31) 



is the distance of atom i from the center of mass Rem of 
the cluster, 



Rr 



El 



(32) 



The force constant was set to kcont = 10 kcal/mol A^^ 
and the distance Rconi ^^ which the atoms start to feel the 
confining potential was set to i?conf = 18 A, sufficiently 



large to accommodate typical shape fluctuations of the 
cluster. Due to this confining potential the cluster is 
effectively in equilibrium with the vapor phase. 

For the calculation of the free energy F{r) as a function 
of the distance of the center of charge from the center of 
mass we have performed a series of calculation in which 
a parabolic biasing potential 



VB(a;) 



-(^-■Rbias)^ 



(33) 



was added to the potential energy of the system. The 
biasing potential depends on the configuration of the sys- 
tem only through 



Rf 



(34) 



the distance of the center of charge from the center of 
mass. Forces resulting from the bias potential given in 
Equ. (j33p can be calculated using the formalism pre- 
sented in Sec. IIIII With a sufficiently large value of the 
force constant fcbias, a biasing potential of this form forces 
the distance r between the center or charge and the cen- 
ter of mass to stay close to i?bias- By combining several 
simulations with appropriately selected values of the pa- 
rameters fcbias and i?bias the complete range of interest 
for the distance r can be sampled. In particular, the 
bias potential coerces the simulation to visit configura- 
tion space regions that would not be visited without bias 
on time scales accessible to a simulation without bias. 
The simulation parameters used in this work are listed 
in Tabic U 

As one is usually interested in the properties of the 
system without the bias, one must correct for its effect 
on the observables 33]. Here we are interested in the 
distribution P{r) of the distance r. 



P(r) = (J[r-r(x)]), 



(35) 



where the angled brackets (• • •) denotes a canonical av- 
erage. In this case, the correction for the bias yields 



P{r) ex cxp[VB{r)/kBT]{S[r - r{x)])u: 



(36) 



where k^ is Boltzmann's constant and (• • •)bias denotes 
an average in the biased ensemble. The distribution 
{S[r — r'(x)])bias of r can be calculated by histogram- 
ming r in the simulation with bias. In principle, the 
complete distribution of r can be determined according 
to Equ. (j36p from one single simulation. In practice. 



however, the distribution P(r) calculated in this manner 
is accurate only in the range of r where in the biased 
simulation sufBcient statistics is accumulated, i.e., in the 
range near i?bias- To determine the distribution P{r) 
over a wider range of distances r one has to combine dis- 
tributions obtained from various simulations with appro- 
priate values of -Rbias and fcbias- Systematic procedures 
have been developed for this purpose [40]. Here, however, 
it proved sufficient to match the distributions obtained 
from the various biased simulations by hand by multi- 
plication of the individual distributions with appropriate 
constant factors. 



V. RESULTS 

A typical configuration of the protonated cluster ob- 
served in a molecular dynamics simulation carried out at 
T = 280 K without bias is shown in Fig. [1] Visual inspec- 
tion of sequences of cluster configurations does not show 
any indication of crystallinity and indicates that proton 
transfer and the reorganization of the hydrogen bond net- 
work occurs on a picosecond time scale. The qualitative 
perception that the cluster is in its liquid state can be 
made quantitative by considering appropriate time cor- 
relation functions. For the dynamics of hydrogen bonds, 
the time correlation function 



CHB(i) 



{h,,{0)h,,{t)) 



(37) 



measures the conditional probability that a hydrogen 
bond exists between two particular water molecules at 
time t provided it existed at time 0. This time corre- 
lation function was introduced by Luzar and Chandler 
to study the dynamics of hydrogen bonds in bulk water 
|4ll l4a | . In the above expression the indicator function 
hij is unity if there is a hydrogen bond from molecule i 
to molecule j and it vanishes otherwise. Here, a hydro- 
gen bond is defined to exist if the Rqo distance is less 
than 3.5 A and the HOO-angle is less than 30°. The be- 
havior of Chb (i) shown in Fig. [2] for the cluster is very 
similar to that of bulk liquid water at room temperature 
|4l| indicating that the kinetics of hydrogen bonds in our 
cluster occurs at approximately the same time scales as 
that in the liquid phase. 

An analogous time correlation function can be used to 
quantify the kinetics of proton transfer: 



Cp(i) = 



{h{Q)K{t)) 



(38) 



The indicator function hi (t) is unity if oxygen i is the hy- 
dronium oxygen in the EVB-state with the largest weight 
at time t and it vanishes otherwise. The correlation func- 
tion Cp(i) measures the conditional probability that the 
excess proton localized near oxygen i at time is still 
there at a time t later. Similar time correlation functions 
have been previously used by Vuilleumier and Borgis [17[ 




FIG. 1: Cluster configuration from a molecular dynamics run 
at T = 280 K without bias. The hydronium oxygen of the 
EVB-state with the largest weight is depicted in yellow and 
hydrogen bonds are shown as blue dashed lines. Typically, the 
excess proton exists as an Eigen-cation at the cluster surface, 
but other configurations also occur. 



and Schmitt and Voth Tl| to study proton transfer in 
bulk water. The form of the correlation function Cp(t) 
results from an interplay between the hopping of the ex- 
cess charge from one molecule to another and the diffu- 
sion of water molecules in the cluster. The similarity of 
the two time correlation functions for hydrogen bonding 
and proton transfer shown in Fig. [5] is a reflection of 
the crucial role played by hydrogen bonds in the proton 
transfer process |43j . 
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FIG. 2: Time correlation function C{t) = {h{0)h{t))/{h) 
for proton transfer (sohd line) and hydrogen bond dynam- 
ics (dashed line). The horizontal dotted line denotes the 
long-time limit of the correlation function for proton trans- 
fer, C — > 1/Nw- The inset shows the same curves with a 
logarithmic scale on the y-axis. 

As exemplified by the configuration shown in Fig. [U 
the excess proton is predominantly located at the cluster 
surface. The stabilization of this position is sufficiently 



strong to prevent the excess charge from visiting the clus- 
ter interior even in molecular dynamics simulations that 
are several nanoseconds long. This observation is con- 
firmed by free energy calculations carried out with bi- 
ased simulations as described in Sees. IIIII and IIVI The 
probability distribution Pproton(?') of the distance r of 
the center of charge from the center of mass, calculated 
from several biased simulations according to Equ. p6p . 
is shown in Fig. [3] along with the distribution Pwatcr('') 
of the distance r of water oxygens from the center of 
mass obtained from an unbiased simulation. As can be 
inferred from the figure, the distribution for the proton 
is peaked at a slightly larger distance than that for the 
water molecules and falls off more rapidly particularly for 
smaller distances. 
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FIG. 3; Probability distribution P(r) of the distance r of the 
center of charge (solid line) and of water oxygens (dashed line) 
from the center of mass of the cluster. 



The form of the probability distributions shown in 
Fig. [3] is strongly influenced by the configuration space 
available for different values of the radius r. To sep- 
arate this purely geometric factor from less trivial ef- 
fects we have calculated normalized densities n{r) ob- 
tained by comparison with the corresponding distribu- 
tions for uniform systems. In particular, the normalized 
density nwatcr('') for water molecules was obtained by 
normalization with the distribution expected in a uni- 
form system with the number density pbuik — 55.5mol/l 
of bulk water, nwatcr(?') = A^Pwator(?')/47rr2pbuik- In 
other words, nwater('') is the observed number density 
at distance r from the center of mass compared to the 
number density of the bulk. In the core of the clus- 
ter, i.e., for r < 7 — 8A, the water density is approxi- 
mately constant and slightly exceeds the density of bulk 
water. Between 8 A and 12 A the water density de- 
creases smoothly to zero. This water density profile indi- 
cates that while the cluster fluctuates in shape it remains 
mostly compact. For the center of charge the normalized 
density nproton(?') = Pproton{r) / 'i-rrr'^ Puniiorm is obtained 
by comparison with a uniform proton density Puniform 
corresponding to one single proton in a sphere of radius 
i?o = 10.8 A, which contains 99% of the proton density. 
The proton density peaks at about 9 A, approximately 



where the water density begins to decrease from its bulk 
value. In the cluster core, out to a distance of about 5 
A, the proton density is essentially constant at a value 
that is smaller than the peak density by a factor of about 
40. Note that even though this factor is not exceedingly 
large, spontaneous excursions of the proton to the clus- 
ter center are very rare. One can estimate that during 
a bias-free MD-simulation with a duration of 1 ^s the 
protonic defect would spend less than 50 ps in the core 
region closer than 1 A to the center of mass of the clus- 
ter. With our formalism, statistically adequate sampling 
of all proton positions including the cluster core can be 
achieved using much shorter simulations. 
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FIG. 4: Normalized density n{r) of the center of charge (solid 
line) and of water oxygens (dashed) as a function of the dis- 
tance r to the center of mass. 

From the normalized proton density npi-oton(?') the po- 
tential of mean force 



W{r) = -fcBrln71proton(?'), 



(39) 



follows. The potential of mean force, which is the re- 
versible work required to change the distance r of the 
proton defect from the center of mass of the cluster, dif- 
fers from the free energy F(r) — —k^T In Ppj-otonif) by 
an entropic term that depends on the surface area of a 
spherical shell with radius r. The potential of mean force 
for our 128-molecule cluster is depicted in Fig. [51 where 
W{r) was normalized such that it vanishes at its mini- 
mum. To carry the proton from the surface of the cluster 
at r « 8 A into its interior a reversible work of about 4 
kBT must be expended. Within a core with a radius 
of about 4 A the proton can be translated at no cost. 
This stabilization of the surface position is very similar 
in magnitude to the stabilization of an excess proton near 
a carbon nanotube membrane with respect to the bulk 
found in recent simulations [25]. The preferential loca- 
tion of the excess charge near surfaces was noted in earlier 
simulations _44] and is consistent with experimental data 
[4a | . From the viewpoint of continuum electrostatics this 
is unexpected because a charge buried deep in the droplet 
is solvated in a more favorable way. The molecular de- 
tails of the solvation structure of the excess proton, how- 
ever, lead to a preferred surface position, presumably due 



to the strain exerted on the surrounding hydrogen bond 
network when the excess proton is in the cluster core. 
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FIG. 5: Potential of mean force W{r) — —/cBTln n(r) as a 
function of the distance r of the center of charge from the 
center of mass of the cluster (solid line). The potential of 
mean force is shifted such that its minimum value is zero. The 
dotted line indicates the potential of mean force obtained from 
an incorrect biased MD-simulation, in which the first term on 
the right hand side of Equ. (|8]) was neglected (see main text). 

One may be tempted to simplify the calculation of bias 
forces by neglecting those force contributions that arise 
from the derivatives dcoi/dx°' of the ground state expan- 
sion coefficients cq; on the right hand side of Equ. ([8]). 
In this approximation, the atomic forces resulting from 
the bias arc simply obtained as weighted average with 
weights Cgj of the bias forces in the individual diabatic 
states. We have inspected the magnitudes of these two 
force terms for the bias potentials used in this study and 
we have found that both forces are of similar magnitude. 
Neglecting the derivatives dcoi/dx" is therefore unjusti- 
fied. To exemplify the effect of this approximation we 
have calculated the potential of mean force W{r) from 
an incorrect biased simulation with fcbias ~ 2.0 kcal/mol 
A~^ and i?bias = 0.0 A, in which the first term on the 
right hand side of Equ. ([5]) was neglected (of course, in 
such a calculation energy is not conserved even without 
thermostat) . The potential of mean force resulting from 
this calculation, shown as dotted line in Fig. [5l strongly 
deviates from the correct W{r). A complete force calcu- 
lation using the formalism of Sec. Illll is therefore crucial 
for a correct free energy calculation. 

It is interesting to compare the solvation structure of 
the excess proton at the surface and in the core of the 
cluster. A typical configuration from a biased simulation 
in which the excess proton is near the center of mass of 
the cluster is shown in Fig. [6] Comparison of configu- 
rations with the proton at the surface and in the core 
reveals that the local solvation structure of the excess 
proton is similar in these two cases. In both cases the 
excess proton is preferentially found in the Eigen struc- 
ture, observed also in bulk water f8|, in which a central 
hydronium ion donates strong hydrogen bonds to three 
surrounding water molecules. 




FIG. 6: Cluster configuration from a molecular dynamics sim- 
ulation in which the center of charge was forced to stay close 
to the center of mass by a bias potential with -Rbias = 0.0 A 
and fcbias = 10.0 kcal/mol A~^. As in Fig. [T]the hydronium 
oxygen of the EVB-state with the largest weight is depicted in 
yellow. Also in the interior of the cluster the protonic defects 
exists preferentially as Eigen cations. 



This qualitative picture is confirmed by a calculation 
of the distribution of the proton coordinate Ar. For the 
calculation of this coordinate one first determines the hy- 
dronium oxygen Oi in the state with the largest weight 
cf. Oxygen O2 is then the oxygen atoms closest to oxy- 
gen Oi. The proton coordinate Ar is calculated from the 
distances of Oi and O2 to the hydrogen H participating 
in the hydrogen bond from Oi to O2: Ar = rH02 ~ ''HOi • 
Distributions of Ar for configurations with the proton at 
the surface and in the core are shown in Fig. [T] The dis- 
tribution of Ar for the core region was calculated from a 
biased simulation using a bias potential with i?bias — 0.0 
A and /cbias = 10.0 kcal/mol A^^. In this biased sim- 
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FIG. 7: Symmetrized distribution for the proton coordinate 
Ar for the unbiased cluster (solid line) and the cluster in 
which the center of charge was biased to stay close to the 
center of mass to the cluster (dashed line). Also shown is the 
distribution of Ar in bulk liquid water at T = 300 K (dotted 
line) gi]. 



ulation the center of charge fluctuated around the cen- 
ter of mass with a width of (?'^)bias — 0.1689 A^. As 
can be inferred from the figure, both asymmetric Eigen- 
hke configurations with |Ar| « 0.35 A and symmetric 
Zundcl-hke configurations with |Ar| « 0.0 A occur and 
in both cases, the Eigen-Hke configurations are preferred. 
In the core of the cluster, Eigen-like configurations occur 
slightly more frequently than at the cluster surface. The 
distribution of Ar in the core is very similar to that ob- 
served in liquid water at a slightly higher temperature of 
T — 300 K shown in Fig. [7] as a dotted line. Thus, as far 
as the local structure of the excess proton is concerned, 
the cluster core provides essentially the same solvating 
environment as the bulk liquid. 



VI. CONCLUSIONS 

In the present paper we have presented a new algo- 
rithm for the calculation of forces in EVB models where 
a bias potential is applied to the center of charge, ef- 
fectively the location of the excess proton. Such bias 
potentials are necessary if one is interested in exploring 
unlikely but important configurations that occur, for in- 
stance, during transfer events. The procedure, based on 
first order perturbation theory, is easy to use and does not 
require any significant additional computational effort as 
all needed quantities are already available from the reg- 
ular EVB force calculation. Bias forces calculated in this 
way for the are rigorously correct (for the EVB-model) 
because they have been derived using perturbation the- 
ory rather than by finite differences. 

The algorithm presented in this paper can be used to 
study the free energetics of proton transfer quantitatively 
in all systems amenable to an empirical valence bond de- 
scription including proton transfer in biological systems 
and acid-base chemistry. As an illustrative example, we 
have calculated the reversible work required to carry an 



excess proton from the surface of a small water droplet 
to its center. The excess charge is preferentially located 
at the surface where it is stabilized by about 4 fceT' with 
respect to a position in the core of the cluster. This pro- 
clivity of the proton to occur near interfaces, which has 
been previously observed in clusters [23, |23| , at planar 
liquid- vapor interfaces [44], and near hydrophobic mem- 
branes [25]] , is unexpected from continuum electrostatics 
and is most likely due to the strain exerted by the ex- 
cess charge on the surrounding hydrogen bond network. 
Analogous behavior has also been observed for other ions 
in clusters [13, IH, Si] . 

Another important application of the formalism devel- 
oped in this paper is the calculation of rate constants 
for proton translocation for instance through membranes 
[2g | using the reactive flux formalism of Bennett and 
Chandler [50, |51| , in which two separate simulations are 
carried out. In both parts the formalism of this paper 
can be applied. In the first step, one determines the 
reversible work required to move the proton across the 
membrane. This can be done efficiently with a molecular 
dynamics simulation with a bias on the center of charge. 
In a subsequent step, the so called transmission coeffi- 
cient is determined by releasing dynamical trajectories 
from a dividing surface located at the free energy maxi- 
mum. The initial conditions needed for such a procedure 
can be generated efficiently by constraining the center of 
charge to be located on the dividing surface by apply- 
ing an appropriate bias potential acting on the center of 
charge. 
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